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SOLUTION  OF  THE  BOLTZMANN  EQUATION  USING 
AN  ENERGY  GROUP  METHOD 

I.  Introduction 

In  certain  regimes  of  beam  propagation,  it  has  been  shown  that  high 
energy  plasma  electrons  whose  mean  free  path  is  long  (but  not  infinite)  may 
play  an  important  role.  Such  electrons  may  be  delta  rays  (beam-generated 
secondaries)  or  thermal  electrons  accelerated  by  large  electric  fields. 
Since  high  energy  electrons  produce  nonlocal  transport,  they  are  not 
correctly  represented  by  eitner  Ohm's  conductivity  law  or  by  fluid 
equations.  The  correct  representation  is  the  linear  Boltzmann  equation  - 
linear  because  electron-electron  collisions  may  be  neglected  but  collisions 
of  electrons  with  nearly-stationary  ions  must  be  included.  Needless  to 
say,  the  linear  Boltzmann  equation  is  also  the  central  equation  in  many 
other  scientific  disciplines.  However,  the  Boltzmann  equation  is  a  six¬ 
dimensional  integro-differential  equation  in  phase-space  (plus  time)  and 
can  only  be  solved  by  approximate  methods.  Many  attempts  have  been  made 
at  solving  that  equation  and  various  approximations  have  been  made.  Within 
the  beam  propagation  context,  analysis  of  nonlocal  transport  by  high  energy 

plasma  electrons  has  depended  primarily  a  'onte  Carlo^  and  particle 

2 

simulation  techniques,  which  are  strongly  unstrained  by  computer  time 

3 

and/or  storage;  or  fluid  approximations,  which  are  very  approximate;  or 
spherical  harmonic  expansions  of  the  Boltzmann  equation. The  latter 
technique  has  been  highly  successful  in  deducing  atomic  and  molecular 
properties  from  swarm  experiments,  and  it  has  been  argued  that  even  two- 
term  spherical  harmonic  expansions  are  valid  even  for  highly  anisotropic 
cases. ^  Indeed  there  are  problems  and  particular  cross-sections  for  which 
two-term  expansions  are  accurate,  but,  in  general,  they  are  very 
inaccurate,  even  divergent,  in  the  anisotropic  limit  of  nearly 
collisionless  streaming. 

Manuscript  approved  August  9,  1989. 
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In  this  report,  ve  describe  and  explore  a  new  technique  which  is 
based  on  choosing  the  energy  of  particles  as  an  independent  variable.  This 
selection  is  suggested  by  the  strong  dependence  of  high  energy  electron 
mean  free  paths  and/or  scattering  kernels  on  the  electron  speed.  Before 
applying  this  method  to  full  scale  problems,  ve  have  applied  it  to  various 
collisionless  test  problems  ranging  from  the  simplest  (force-free,  free- 
expanding  electrons),  to  more  complicated  ones  like  shear  flow  in  one 
dimension.  These  test  problems  have  two  characteristics: 

-  They  allow  a  comparison  with  exact  analytical  solutions 
which  can  be  obtained  directly  from  the  original  linear  Boltzmann  equation. 

-  Because  these  tests  are  collisionless,  they  constitute 
severe  tests  of  the  numerical  solution  since  no  physical  smoothing  of  the 
solution  is  provided  in  the  equations. 

The  outline  of  this  report  is  the  following.  In  the  first  section, 
we  illustrate  problems  with  some  of  the  previous  methods  used  to  solve  the 
linear  Boltzmann  equation.  In  the  second  section,  we  present  the  equations 
and  the  approximation  we  make  to  solve  these  equations.  In  the  third  part, 
the  numerical  scheme  is  presented.  In  the  ensuing  sections,  we  report 
successively  on  the  following  collisionless  test  problems:  free-streaming, 
free-streaming  in  the  direction  parallel  to  the  flow  and  perpendicular  to 
the  flow  and  finally  shear  flow.  A  review  of  these  results  and 
recommendations  for  future  work  will  then  follow  in  the  final  section. 

II.  Example  of  Problems  Associated  with  Approximations  in  Solving 

Boltzmann  Equation 

In  this  section,  we  consider  a  simple  problem,  i.e.,  collisionless 
steady-state,  one-dimensional  cold  flow  of  electrons  subject  to  a  constant 
electric  field  E  .  In  this  case,  the  Boltzmann  equation  reduces  to 
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with  the  boundary  conditions 

f(z  =  0,  v)  =  5(vx)6(vy)5(vz  -  vq), 


(2a) 


f(z  -»  ®,v)  =  0  for  <  0. 


(2b) 


We  assume  that  qE  >  0,  so  that  electrons  are  accelerated  toward  z  =  The 
exact  solution  is 


f(z,v)  =  S(vx)S(Vy)S 
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where  9  is  the  step  function.  Taking  moments,  we  find  that  the  density 

N(z)  and  fluid  velocity  U  (z)  are  given  by  the  first  column  of  Table  1, 

z 

2  2 

wherein  V  is  defined  by  V  =  U  +  2qEz/m. 
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We  have  also  worked  out  the  values  of  these  macroscopic  quantities  by 
solving  the  two-polynomial  and  four-polynomial  expansions  of  the  Boltzmann 
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equation.  It  is  seen  in  Table  1  that  NU2  is  given  exactly  by  the 

polynomial  expansions,  but  that  with  increasing  z,  the  tvo-polynomial 

expansion  shows  N  increasing  monotonically  to  infinity  whereas  the  exact 

solution  has  N  decreasing  monotonically  to  zero,  and  conversely  for  U  . 

Thus,  the  polynomial  solution  could  hardly  be  more  wrong.  Even  more 

disturbing,  the  four  polynomial  expansion  worsens  the  error. 

(Incidentally,  if  an  initial  value  spatially-uniform  problem  were  solved, 

rather  than  a  steady  state  problem,  N  would  be  given  correctly  by  the 

polynomial  approximations  but  U  and  NU  would  be  totally  incorrect.)  The 

z  z 

error  can  be  traced  to  the  fact  that  the  few-polynomial  expansions  assume 
near-isotropy,  and,  therefore,  incorrectly  heat  the  transverse  dimensions. 
Equally  disturbing  is  the  absence  of  nonlinear  conductive  terms  in  the 
polynomial  formulation. 

The  method  we  propose  here  is  an  alternative  type  of  expansion  which 
also  focuses  on  density  and  current,  the  quantities  of  principal  interest, 
and  also  retains  energy  as  an  independent  variable,  but  which  becomes  ex^ct 
in  both  the  near  isotropic  limit  and  the  cold  fluid  nearly  collisionless 
limit. 

III.  Equations 

The  method  starts  with  the  linear  Boltzmann  equation 

It-  +  -  ‘  lx  +  m  [qE(x,  t )  +  qv  x  B(x,  t)  j  •  =  Jd3v'K(v'  ,  v)  f  ( x'  ,  v'  ,  t ) , 

(4) 

where  f  is  the  distribution  function  in  space  and  velocity  and  K  is  the 
scattering  kernel.  Now,  by  taking  moments  over  velocity  angles  while 
keeping  the  speed  w,  i.e.,  the  kinetic  energy,  as  an  independent  variable, 
the  Boltzmann  equation  yields  two  coupled  equations: 
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(5b) 

where  n(x,t,v)  is  the  "microscopic"  density  of  particles  with  speed  v,  and 
u  is  the  mean  velocity  of  the  class  of  electrons  with  speed  w,  defined 
respectively  as 

n(x,t,w)  =  Jd3v  5  ( | v |  -  w)  f  (x,v,t),  (6) 

n(x, t , w)  u  (x,t,w)  =  jd3v  S  (|v|  -  w)  v  f(x,v,t)  (7) 

(w,w')  and  (w,w' )  are  integrated  scattering  kernels  and  are  defined 
as 


-a 

Kq  (w,w')  =  2n  J  ax  sin  X  K  (w,v',X) 


(8) 


.11 

K1  <v,v')  =  2 It  Jd X  sin  X  cos  X  K  (v,w', 


X) 


(9) 


where  X  is  the  angle  between  v  and  v' . 

To  close  the  hierarchy  of  moment  equations,  one  approximation  is 
made:  we  assume  that  the  microscopic  pressure  tensor  (for  particles  with 
speed  w) 


£(x,t,w)  =  m  Jd3v  5  (|vj  -  w)  £v  -  u(v)j  £v  -  u(w) 


f(v)  (10) 
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can  be  replaced  by  an  isotropic  scalar  value 

p(x,t,w)  =  j  n(v)  [v-  -  u2(w)j  (11) 

which  is  then  uniquely  specified  by  the  requirement  that  the  particles  in 
question  have  speed  w. 

Let  us  note  that  this  form  of  the  pressure  represents  exactly  the  sum 
of  the  three  diagonal  components  of  p(z,t,w)  i.e.,  it  does  not  mix 
streaming  energy  with  thermal  energy  (as  do  low-order  spherical  harmonic 
expansions  of  the  Boltzmann  equation,  for  example).  This  approximation 
would  appear  to  be  a  relatively  weak,  one  compared  to  those  that  are  made  in 
deriving  macroscopic  fluid  approximations;  however,  we  shall  see  that  it 
does  lead  to  some  serious  problems.  This  technique  reduces  the  six¬ 
dimensional  Boltzmann  equation  to  a  pair  of  four-dimensional  equations  and 
a  great  simplification.  Ve  shall  concentrate  in  this  report,  however,  on 
potential  problems  and  deficiencies,  and  how  we  have  explored  them.  In 
particular,  the  isotropic-microscopic-pressure  assumption  is  weakest  in  the 
completely  collisionless  limit.  In  this  limit,  it  is  often  possible  to 
write  down  exact  solutions  of  Eq.  (4)  in  closed  form  in  terms  of  simple 
integrals  over  the  deterministic  orbits.  Ve  have  thus  been  able  to  compare 
these  exact  solutions  to  those  of  our  model.  Let  us  keep  in  mind  that  we 
are  discussing  a  model  for  the  Boltzmann  equation,  rather  than  a  precision 
approximation  scheme.  (Even  so,  this  would  be  an  improvement  in  the 
present  state-of-the-art.)  It  is  important  to  note  that  our  real  interest 
is  in  calculating  the  macroscopic  electron  density  and  current,  for  use  as 
sources  in  the  field  equations:  getting  these  quantities  right  should  be 
easier  than  calculating  all  the  microscopic  distributions  accurately. 
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Displaced  Maxwellian 

We  may  first  ask.  hov  veil  Eq.  (11)  represents  p(z,t,v)  in  simple  test 
cases,  like  displaced  Maxwellian  velocity  distributions, 


f(v)  =  exp 


-( v-U)  /v 


2‘ 

thj  ‘ 


(12) 


In  Fig.  1,  we  show  the  exact  values  of  p.  and  p. ,  obtained  by  substituting 
Eq.  (12)  in  Eq.  (10),  and  compare  them  with  the  assumed  isotropic  value  in 
Eq.  (11).  We  note  that:  (a)  p|(  and  p^  are  correct  to  first  order  in  u/w 
in  the  isotropic  limit  u  -»  0.  (b)  •*  0  and  p,  ->  0  in  the  streaming  limit 

u  -»  v.  Here  p  -*  0  as  well.  However,  the  exact  solution  ;s  not  isocropic 
in  this  limit;  rather  p^  /p^  •*  0.  We  argue  that  our  representation  is 

adequate,  in  that  plays  a  vanishingly  small  role  in  this  limit,  where  the 

"  2 
dynamics  are  dominated  by  the  streaming  energy  1/2  mU  .  Thus,  our 

approximation  is  valid  to  about  20%,  where  it  matters. 


Anomalous  Cases 

It  is  easy  to  conceive  of  particular  velocity  dis tributions  for  which 
isotropic-microscopic-pressure  is  not  a  good  approximation.  We  have 
already  seen  one,  but  in  a  limit  where  the  value  of  the  pressure  is  not 
important.  Consider  next  a  counterstreaming  case, 

r(v)  =  a  S(v-U)  +  b  o( v-U) .  (13) 

2  2 

Here  p^  =0  and  therefore  p^  (v)  =»  n  ( v4"  -  u“).  On  the  other  hand, 
consider  a  distribution  which  is  initially  spatially  localized  eo  a  thin 

sheet  at 

z  =  0  and  spreads  out  through  collisionless  free-streaming, 

f(z,v,t  =  0)  =  S(z)  exp(-  v2/v^h),  (I4) 
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(15) 


which  givj - 

f  ( 2  ,  v ,  )  =  (r  -  vz  t)  exp(-  v2/v2h). 

Since  v^  is  here  uniquely  specified  as  a  function  of  2,  we  have  =  0  and 
therefore , 

p^(w)  =  1/2  n  (v2  -  u2).  (16) 

These  examples  illustrate  the  point  that  there  is  no  unique 
microscopic  pressure  which  can  correctly  represent  all  the  structure  that 
can  evolve  in  six-dimensional  phase  space,  especially  in  the  coliisionless 
limit.  Indeed,  other  macroscopic  models  inevitably  suffer  from  the  same 
deficiency:  for  example,  the  spherical  harmonic  expansion  of  the  Boltzmann 

equation  would  require  an  enormous  number  of  polynomials  to  represent 
either  of  the  cases  described,  and  this  would  amount  essentially  to  fully 
resolving  3-D  velocity  space.  We  do  not  regard  this  as  a  fatal  flaw  in 
our  method,  inasmuch  as  one  of  our  objectives  is  modest:  to  calculate  the 
macroscopic  density  and  current  to  reasonable  qualitative  accuracy  (10- 
302).  To  this  end,  our  assumption  has  the  advantages  of  simplicity  and 
centrality,  i.a.,  it  appears  equally  likely  to  err  in  either  direction.  We 
have  abandoned  earlier  attempts  to  define  an  anisotropic  pressure 
prescription  that  is  more  accurate  than  Eq.  (11)  in  (for  example)  the 
single-streaming  limit.  Such  formulas  are  inevitably  found  to  be  less 
accurate  in  other  limits,  and  also  sacrifice  mathematical  simplicity. 
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Sound  Waves 


A  collisionless  field-free  gas  supports  no  sound  waves.  However, 
Eqs.  (5a),  (5b)  and  (11)  do  support  a  sound  wave,  with  speed  cs  (v)  = 

7  2  1/2 

(w“  -u  )  ,  within  each  energy  population.  When  v  is  integrated  out 

obtain  macroscopic  quantities,  these  waves  phase  mix  and  damp  away  as 
2  7  2 

exp(-  k  v^  t  /12)  (for  the  case  of  a  Maxvell-Bol tzmann  distribution); 
furthermore,  the  waves  are  damped  if  there  are  collisions  or  fields 
coupling  different  energy  groups. 


to 


IV.  Numerical  Solution 

In  the  case  of  one  spatial  dimension,  the  code  is  similar  in  structure 
to  a  2-D  time-dependent  magnetohydrodynamic  code  and  consists  of  solving 
Eqs.  (5).  There  is  no  energy  equation  since  the  pressure  is  specified  in 
closed  form.  In  a  first  stage,  the  code  has  been  implemented  in  the 
simplest  possible  geometry,  i.e.,  cartesian  1-D  in  space  but  is  2-D-like 
because  w  (i.e.,  energy)  is  retained  as  an  independent  variable.  These  two 
dimensions  are  not  treated  similarly  numerically;  for  example,  the  gridding 
is  linear  in  space  and  logarithmic  in  energy.  In  collisionless  cases  and 
in  the  absence  of  electric  fields  as  in  the  test  cases  considered  in  this 
report,  there  is  no  exchange  between  the  different  energy  groups  and  each 
energy  group  is  completely  independent  of  its  neighbors.  Operator¬ 
splitting  has  been  used  to  treat  the  two  dimensions  and  since  the 
microscopic  densities  convect  in  the  same  way  as  fluid  densities,  Flux- 
Corrected  Transport"7  has  been  chosen  for  the  convection  algorithm.  Also, 
the  equations  have  been  recast  in  conservation  form  and  we  note  that  the 
equivalent  convection  velocity  in  the  microscopic  velocity  equation  is 
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where  the  symbols  have  been  defined  previously  and  where  y  is  the 
relativistic  factor. 

The  important  output  of  the  code  is  not  the  "microscopic"  (i.e., 
w-dependent)  quantities  which  are  carried  from  time  step  to  time  step  but 
rather  the  macroscopic  variables  such  as  the  macroscopic  density  N(z,t) 
and  the  average  macroscopic  velocity  V(z,t).  (We  use  capital  letters  to 
denote  these  macroscopic  quantities.)  These  are  integrals  over  w  of  the 
microscopic  variables  n(z,t,w)  and  u(z,t,w).  As  mentioned  before,  it  is 
those  macroscopic  variables  which  may,  in  turn,  be  used  as  input  into  other 
codes  or  which  can  be  compared  to  other  fluid-type  calculations. 

A  major  difference  between  this  code  and  a  fluid  code  appears  in  the 
time  step  constraint.  In  an  explicit  hydrodynamic  code,  the  Courant  time 
step  condition  can  be  written 

St  <  0.5  j-  p;--  ,  (18) 


where  v  is  the  fluid  velocity  and  cs  is  the  speed  of  sound.  As  seen 
before,  the  sound  speed,  as  derived  from  the  isotropic-microscopic- 
pressure,  is  equal  to 


1/2 


afp(y  --U  )] 

3p 


.  2  2. 

=  (V  -u  ) 


1/2 


(19) 


We  see  tnat  speed  of  sound  is  a  function  of  the  energy  group  and  has  a 
maximum  for  each  group.  So,  as  the  energy  of  the  group  increases , the  time 
step  decreases.  If  there  is  a  need  to  include  very  high  energy  groups  in 
a  calculation,  the  time  step  may  become  very  small  Since  the  higher 
energy  groups  contain  fewer  particles,  we  would  be  in  a  situation  where  the 
regions  of  smaller  densities  would  control  the  time  step.  This  situation 
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is  clearly  undesirable  numerically  and  moreover,  this  time  step 
restriction  is  not  physical  in  origin  since,  for  a  collisionless  fluid, 
there  is  no  true  sound  speed.  The  situation  arises  from  the  approximate 
(isotropic)  expression  used  for  the  microscopic  pressure. 

Typically,  there  are  many  fever  high-energy  particles,  and  these 
particles  have  longer  mean  free  paths  and  thus  have  a  coarser  spatial 
structure.  Therefore,  a  variable  space  gridding  has  been  devised  and 
implemented  for  these  groups  of  particles  to  speed  up  the  calculations. 

The  principle  behind  this  variable  gridding  is  the  following:  above  a 
predetermined  velocity  not  to  exceed  twice  the  thermal  velocity,  the 
number  of  points  in  the  spatial  dimension  is  divided  by  two.  For  groups 
with  speed  greater  than  2  u^,  the  number  of  grid  points  is  divided  by  two 
once  more  and  so  on.  In  this  fashion,  the  overall  time  step  criterion  is 
applied  to  the  velocity  uc  which  is  of  the  order  of  the  thermal  velocity 
and  which  is  much  smaller  than  the  largest  velocity  in  the  calculation. 

The  total  number  of  grid  points  in  the  lowest  energy  group  does  not 
necessarily  have  to  be  a  power  of  two  since  the  last  spatial  point  can  be 
dropped  out  of  the  next  regridding.  There  is  a  limit  to  the  number  of 
divisions  by  two  which  can  be  accomplished  practically  since  we  must  be 
left  in  the  end  with  at  least  half  a  dozen  grid  points  or  so  for  the 
highest  energy  group.  Because  of  that  consideration,  the  time  step  may 
not  be  controlled  by  the  lowest  energy  group  but  still  by  the  highest  one 
on  a  much  coarser  grid,  however. 

Several  comments  can  be  made  about  this  approach.  This  technique  can 
easily  be  extended  to  multi-spatial  dimensions  and  in  tne  case  of  two 
spatial  dimensions,  for  example,  the  cell  becomes  4  times  larger  (a  factor 
of  two  in  each  direction)  and  in  3-D,  the  cell  becomes  8  times  larger  at 
each  regridding.  This  technique  also  requires  a  repacking  of  the  arrays  in 
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energy  space.  This  repacking  can  be  done  once  at  the  beginning  of  the 
calculation  and  has  to  be  unfolded  at  each  reconstruction  of  the 
macroscopic  variables,  namely  for  every  time  step  when  output  is  desired. 
That  unfolding  requires  interpolations.  Also,  if  there  are  particle 
exchanges  between  energy  groups  (when  collisions  or  electric  fields  are 
present),  care  must  be  applied  since  spatial  resolution  for  adjacent 
energy  groups  may  differ  by  a  factor  of  two. 

Finally,  instead  of  increasing  the  spatial  grid  for  the  higher  energy 
groups,  the  time  evolution  of  these  groups  could  have  been  computed  on  a 
shorter  time  scale,  basically  subcycling  the  calculation  for  these  groups. 
It  is  easy  to  see  that  such  a  procedure  would  involve  many  more 
calculations  than  the  approach  described  above. 


V.  Test  Problems 

1.  Force-free  free-streaming  Maxwellian  distribution 
In  a  first  series  of  tests,  we  have  looked  at  the  collisionless  case 
in  which  the  electron  distribution  is  initialized  as  a  thin  (but  not 
infinitely  thin)  sheet  in  space  with  no  initial  mean  velocity  around  a 
Maxwellian  velocity  distribution.  The  initial  distribution  function  is 
simply  expressed  as 

2  2 

f(z,v,0)  =  exp  f-  exp  [-  (20) 

*  4 


The  solution  for  f  at  all  times  is  given  by: 


f(z,v, t) 


(21) 
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To  implement  the  present  model  for  any  given  problem,  the  initial 
microscopic  density  and  velocity  n(z,0,v)  and  u(z,0,w)  have  to  be 
calculated  from  f(z,v,0).  The  analytical  solutions  for  n(z,t,v)  and 
u(z,t,v)  are  then  calculated  by  taking  integrals  over  v  and  compared  to  the 
solutions  of  the  model  Eqs.  (5).  In  fact,  these  steps  represent  a 
significant  amount  of  algebra  and  since  these  manipulations  do  not 
constitute  the  main  goals  of  this  study,  they  will  not  be  presented  in 
detail  here.  For  the  very  simple  example  under  consideration,  even  the 
pressure,  as  defined  in  Eq.  11,  can  be  computed  exactly  in  terms  of  a 
simple  integral: 


Figures  2  and  3  shov  n(z,t,v)  and  u(z,t,w)  for  two  values  of  w 
(0.6  v  and  1.6  v  at  t  =  2.96  ns  for  6  =  10  dx  where  dx  is  the  spatial 
grid  step.  For  most  of  the  test  problems,  we  have  chosen  dx  =  0.1  cm, 

3 

v  ^  =  10  cm/s,  the  total  length  of  the  grid  of  the  order  of  2  cm  (thus,  a 
perturbation  at  mid-grid  will  then  take  on  the  order  of  10  ns  to  travel  all 
the  way  to  the  grid  boundary)  and  depending  on  the  problem,  v  ranging  from 
about  0.25  v  ^  up  to  3  or  4  v  Figure  3  shows  a  dip  in  n  near  its 
maximum  value  for  the  larger  energy  group.  This  dip  increases  when  A  is 
reduced  and  for  example,  Fig.  4  shows  n(z,t,w)  and  u(z,t,w)  for  w  =  1.6  v  ^ 
and  for  A  =  3  dx  at  t  =  2.82  ns. 
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Because  of  this  unexpected  feature,  an  exact  analytic  solution  of  the 
model  equations,  i.e.  those  including  the  assumption  of  an  isotropic- 
microscopic  pressure,  was  sought  for  that  simple  case.  Noting  that  each 
w-group  is  independent,  the  equations  can  be  written: 


Bn 

3t  + 


n 


3u 

3z 


0, 


(23) 


3u 

at  + 


1  a  r  (  2 

3n  3z  Lnr 


(24) 


Assuming  that  a  solution  of  the  form 


Twz 

u(z,t)  =  <R( t ) 

z  <  R(t) 

lw 

z  >  R(t) 

(25) 

n(z, t)  =  n.(t) 

z  <  R(t) 

(26) 

exists,  it  can  be  shown  that 


R(t)  =  i  wt,  (27) 

and  finally  the  solution  for  n(z,t)  and  u(z,t)  is  obtained: 


3  3 

t  /t 

0  < 

z  <  jwt 

0 

[iut°  l3 

J  wt 

<  z  <  w  ft 

2  r  1 

b(wt-z)  J 

3  toJ 

1°  Vit  “  I  to)  <  z 

(28) 


u(z, t) 


r 

3  z/t 
w 


1° 


0  <  z  <  j  wt 

J  <  2  <  u  ('  -  J  '„) 
*  ('  - 1 '.)  <  * 

(29) 


This  solution  has  been  plotted  in  Fig.  5  and  shows  a  dip  in  the  middle  for 
each  and  every  group  w!  Thus,  even  for  this  most  simple  problem,  the  model 
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equations  adopted  here  (not  the  numerical  implementation)  give  a  poor 
approximation  to  the  exact  solution  of  the  Boltzmann  equation,  at  least  for 
the  individual  group  v.  The  error  turns  out  to  be  more  than  the  20-30£ 
vhich  vas  thought  of  as  acceptable  at  the  start  of  this  study  and  a  more 
careful  look  at  the  approximation  had  to  be  taken. 

As  indicated  before,  the  only  approximation  vhich  vas  made  in  the 
model  vas  the  assumption  of  an  isotropic-microscopic  pressure.  To  improve 
the  accuracy  for  the  present  test  problem,  ve  have  abandoned  this 
approximation  and  in  its  place  used  a  tensor  pressure  equation.  The 
general  form  for  the  pressure  tensor  equation  is 


3Pij  3  (  )  fjjj.  p  3  Q.  ■ 

3t  +  3x^  lukPijJ  +  Pik  3x^  +  jk  3x^  +  3x^  ^ 


3  Q.  ..  q  E,  3  (1  P .  . i 

3x,  +  m  3v  [v  uk 


To  close  the  moment  hierarchy,  our  next  approach  is  to  drop  the  higher 
order  terms  Q  that  appear  now  in  the  equation.  Further,  using  symmetries 
and  properties  for  the  diagonal  elements,  the  original  model  has  been 
increased  by  tvo  more  equations  (in  each  independent  variable)  and  tvo 
more  quantities  have  to  be  initialized.  These  are  components  of  the 
pressure  tensor  and  more  calculations  have  to  be  performed  before  starting 
a  problem.  However,  once  all  of  these  steps  were  accomplished,  the 
original  problem  vas  attempted  again  for  a  delta- function  in  space 
initially  (the  limiting  case  of  the  problem  vith  vhich  the  original  model 
had  difficulty).  Results  are  shown  for  n(z,t,v)  and  u(z,t,v)  for 
v  =  1.6  v  ^  at  t  =  3.6  ns  in  Fig.  6.  Ve  see  now  that  the  results  are 
quite  acceptable.  The  model  is  not  quite  as  simple  anymore  but  the  results 
are  much  better  than  the  hoped  for  overall  accuracy  of  10-30X. 


2.  Displaced  Maxwellian  distribution  streaming  in  the  direction 
of  motion. 

This  test  problem  is  very  similar  to  the  first  one  and  involves 
mainly  a  change  in  frame  of  reference.  The  v-integrated  macroscopic 
density  and  velocity  N  and  V  should  be  independent  of  the  streaming 
velocity  u  except  for  obvious  displacement  in  z  and  v  .  However,  since  the 
model  is  not  explicitly  frame-independent,  this  test  is  necessary.  The 
initial  distribution  function  for  this  case  can  be  written: 


f(z, v,0) 


2  .  .2 

z  /A 


exp 


(31) 


For  example,  for  this  case,  the  microscopic  density  can  be  shown  to  be  at 
all  times: 
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The  results  for  this  case  depend  on  the  number  of  energy  groups  carried  in 
the  computation  and  it  was  found  that  twelve  groups  were  needed  in  order  to 
reproduce  reasonably  well  the  macroscopic  density  and  velocity.  Figure  7 
shows  the  microscopic  density  and  velocity  at  t  =  3.9  ns  for  the  11 ^ 
energy  group  and  the  macroscopic  density  and  velocity  at  the  same  time. 

Note  that  for  this  problem  t  0  initially  but  P  =  0.  The  agreement 
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for  this  case  is  excellent  between  the  analytical  and  numerical  solutions 
and  this  change  of  reference  frame  due  to  streaming  is  not  a  problem  for 
the  model. 

3.  Displaced  Maxwellian  distribution  streaming  in  the  direction 
perpendicular  to  the  direction  of  motion. 

Analytically,  this  case  is  more  complicated  than  the  previous  one  and 
for  example,  both  components  of  the  pressure  tensor  are  different  from  0 
at  t  =  0.  This  test  problem  was  originally  carried  without  a  pressure 
tensor  formulation.  Actually,  it  was  this  test,  as  well  as  the  first 
test  problem  (force-free,  free-expansion  with  a  localized  distribution  in 
space),  which  showed  enough  discrepancy  with  the  analytical  solution  and 
which  forced  a  re-evaluation  of  the  original  model  as  presented  in  the 
second  section  of  this  report.  Results  for  this  test  are  shown  in  Fig.  8 
at  t  =■  1  ns  both  for  the  macroscopic  and  microscopic  quantities  (for  the 
10th  energy  group).  This  calculation  was  carried  out  with  12  energy  groups 
again  and  for  u  =  2  v  The  agreement  between  numerical  and  analytical 
results  is  not  as  good  as  for  the  previous  case  but  by  comparison  with  the 
scalar  pressure  case  formulation,  it  is  a  much-needed  improvement  and  is 
quite  satisfactory  in  view  of  our  criterion  for  accuracy.  The  numerical 
problem  which  remains  seems  to  be  spreading  of  the  macroscopic  density,  or 
diffusion,  in  the  z-direction  which  increases  with  time.  Numerical  tests, 
e.g.,  changing  the  number  of  energy  groups  did  not  change  the  results  in 
any  significant  way. 

4.  Shear  Flow  Problem 

Increasing  in  complexity,  we  next  look  at  the  collision-free  shear 
flow  problem  characterized  by  the  following  initial  distribution  function 
f: 
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f(z,v,0)  = 
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Its  solution  for  all  times  is  simply 
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and  we  consider  the  special  case  with  the  initial  condition 


uo(z)  = 
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We  can  calculate  the  exact  microscopic  density  and  velocity.  After  some 
algebra,  these  quantities  can  be  shown  to  be: 
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where  9  is  the  Heavyside  function.  The  last  integral  can  be  simplified 
2  2  2 

further  for  z  >  w  t  and  is  written 
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Ue  see  that  n(z,t,w)  is  independent  of  z  and  t  as  it  must  be  and  that  u 
(z,t,v)  is  independent  of  z  and  t  only  when  z  >  vt.  This  is  because,  if  z 
is  positive,  no  particle  with  speed  v  which  originated  in  the  region  z  <  0 
has  yet  arrived  at  z:  thus  there  is  as  yet  no  mixing  of  v  values. 

The  system  of  model  equations  can  be  solved  analytically  for  the 
present  case 


n(z, t , w)  =  cC, 
p32(z,t,w)  =  cC, 
and  can  be  written 
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This  reduces  to  a  wave  equation 
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to  which  the  solution  is  any  function  of  the  form: 


u  (z,t,w)  =  u1  fz  +  -22-  tl  +  u-  (z-  t) 

y'  ’  ’  '  it  ^mn  )  2  t  imn  ) 


(45) 


To  specify  u^  and  Uj,  we  need  initial  values  of  u^  (z,0,w)  and 

[du(z, t,w)/dt]  at  t  a  0.  These  values  can  be  obtained  from  the  analytical 

solution  and  after  some  algebra,  it  becomes 
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where  w(z/t)  is  determined  by  the  condition 
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The  results  for  u  given  by  the  model  are  shown  in  Fig.  9  whereas  the 
analytical  solution  to  the  Boltzmann  equation  is  shown  in  dotted  lines  on 
the  same  figure.  Once  more,  a  discrepancy  exists  between  the  model 
solution  and  the  exact  analytical  solution.  So,  going  back  to  the 
approximations  made  in  the  model,  we  note  that  we  had  closed  the  hierarchy 
by  setting  Q.^  to  0,  on  the  assumption  that  the  Q.^  terms  were  third 
order  in  deviations  from  average  quantities.  But  the  situation  is  really 
more  complicated.  The  difficulty  in  this  shear  flow  problem  originates 


from  the  fact  that  correlations  are  introduced  into  Q...  due  to  the  factor 

1  j  K 
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Ve  have  picked  up  another  tensor  pressure  component  in  the  process,  though 

ve  are  still  solving  a  problem  in  one-dimensional  geometry  only. 

Unfortunately,  ve  are  finding  that  as  the  complexity  of  the  test  problem 

increases,  so  too  does  the  complexity  of  the  model  needed  to  get  a 

satisfactory  solution.  Ve  have  also  found  numerical  problems  in  the  shear 

flow  test  problem,  vhich  also  require  special  treatment.  The  solution  for 

the  shear  flow  problem  for  u  =  0.5  v  ^  is  shown  in  Fig.  10  and  displays 

strong  oscillatory  behavior  around  z  =  0  (where  the  velocity  changes 

its  sign  abruptly).  This  behavior  was  improved  by  dropping  the  FCT 

algorithm  for  P  and  P  ^  component  equations  and  replacing  it  instead  by 
y  z  yy 
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an  upwind  modified  Lax  scheme  for  these  quantities.  In  order  to  update 
these  quantities,  the  operator  9u/3t  had  to  be  replaced  by 

u.n+1  =  0.7  u.n  +  0.15  (u.^  +  u.^),  (52) 

where  n  and  i  are  respectively  the  time  and  spatial  indices.  The 
oscillations  disappeared  until  t  =  3  ns  at  which  time  they  reappeared. 
Another  case  was  run  in  which  the  abrupt  velocity  profile  at  t  =  0  was 
replaced  by  a  smoother  one  of  the  form 

uq(z)  =  uQtanh  (4z).  (53) 

The  value  of  u  was  also  increased  to  1.5  v  ^  and  the  solution  for  that  case 
is  shown  in  Fig.  11  at  t  »  3.9  ns.  The  velocity  profile  steepened  up  in 
time,  trying  to  simulate  the  previous  case.  It  still  showed  some 
oscillatory  behavior  near  z  =  0  at  a  late  time.  Also,  the  solution 
presented  some  small  asymmetry  at  later  times. 

VI.  Conclusion 

In  this  report,  we  have  examined  on  a  variety  of  attempts  to  implement 
a  scheme  to  solve  the  Boltzmann  equation  using  an  energy  group  method.  Our 
first  attempt  to  close  the  hierarchy  of  moment  equations  by  using  a  simple 
approximation  on  the  microscopic  pressure  for  each  energy  group  did  not 
even  hold  for  the  simplest  test  problem,  a  force-free,  free-expanding 
Maxwellian  distribution  function.  Tensor  component  equations  had  to  be 
added  for  that  problem.  Having  done  so,  the  next  test  problems,  involving 
a  free-streaming  distribution  with  mean  velocity  in  the  direction  parallel 
and/or  perpendicular  to  the  spatial  coordinate,  proved  to  yield  a 
satisfactory  solution. 


22 


The  next  problem  which  dealt  with  shear  flow  required  the 
introduction  of  an  extra  equation  for  the  higher  order  heat  flow  tensor. 
Numerical  problems  also  appeared  in  the  treatment  of  a  very  sharp 
transition  in  velocity.  Although  some  of  these  problems  could  have  been 
anticipated  to  some  degree,  our  present  feeling  now  is  that  continuing 
degrees  of  improvement  have  to  be  made  in  order  to  resolve  more  and 
more  complicated  problems.  Our  real  interest  is  in  problems  that  are 
multi-dimensional  spatially,  which  appear  likely  to  introduce  a  good  deal 
more  complexity  into  the  assumption  needed  to  close  the  hierarchy. 

Although  all  the  cases  under  study  were  collisionless,  and  because  of  that, 
were  more  difficult  to  solve  numerically  (since  no  physical  smoothing  was 
present  in  the  equations),  our  confidence  in  the  implementation  of  the 
method  has  not  reached  a  point  where  the  solution  can  be  trusted  in 
untested  physical  situations.  The  method,  which  looked  very  attractive  at 
the  beginning  because  of  its  simplicity,  has  been  shown  to  work  for  the 
problems  studied  after  some  hard  work  and  after  the  addition  of  many 
"corrective"  equations.  The  next  step  of  applying  this  method  to  more 
spatial  dimensions  and  to  more  realistic  problems  has  not  been  pursued  for 
the  moment. 
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n  AT  T  =  2.96  ns  W 


Fig.  2  Microscopic  density  and  velocity  for  a  thin  sheet  (&  =  lOdx) 

force-free  free-streaming  Maxwellian  distribution  function  for 
v  =  0.6  v  ^  energy  group  at  t  =  2.96  ns  using  an  isotropic- 
microscopic  pressure. 
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n  AT  T  =  2.82  ns  W  = 


Fig.  4  Same  as  Fig.  2  for  a  thinner  sheet  (&  =  3  dx)  for  v 


at  t  =  2.82  ns. 


MODEL  SOLUTION 


EXACT  SOLUTION 


89  ns  W 


10  X  107 


Maxwel 


Fig.  9  Model  solution  for  u  for  shear  flow  problem  at  t  =  1.6  ns 


n  AT  T  =  3  94  ns  W 


99  8 
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Fig.  11  Same  di,  Fig.  li)  lui  smoo 
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